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ABSTRACT 

Context. Hot Jupiters exhibit atmospheric temperatures ranging from hundreds to thousands of Kelvin. Because of their large day- 
night temperature differences, condensable species that are stable in the gas phase on the dayside — such as TiO and silicates — may 
condense and gravitationally settle on the nightside. Atmospheric circulation may counterbalance this tendency to gravitationally 
settle. This three dimensional (3D) mixing of chemical species has not previously been studied for hot Jupiters, yet it is crucial to 
assess the existence and distribution of TiO and silicates in the atmospheres of these planets. 

Aims. We investigate the strength of the nightside cold trap in hot-Jupiter atmospheres by investigating the mechanisms and strength 
of the vertical mixing in these stably stratified atmospheres. We apply our model to the particular case of TiO to address the question 
of whether TiO can exist at low pressure in sufficient abundances to produce stratospheric thermal inversions despite the nightside 
cold trap. 

Methods. We perfom 3D global circulation models of HD 209458b including passive tracers that advect with the 3D flow, including 
a source/sink on the nightside to represent condensation and gravitational settling of haze particles. 

Results. We show that global advection patterns produce strong vertical mixing that can keep condensable species lofted as long 
as they are trapped in particles of sizes of a few microns or less on the night side. We show that vertical mixing results not from 
small-scale convection but from the large-scale circulation driven by the day-night heating contrast. Although this vertical mixing is 
not difl'usive in any rigorous sense, a comparison of our results with idealized diffusion models allows a rough estimate of the effective 

vertical eddy difi"usivities in these atmospheres. The value K-- = m /s can be used in ID models of HD 209458b. Moreover, 
our models exhibit strong spatial and temporal variability in the tracer concentration that could result in observable variations during 
either transit or secondary eclipse measurements. Finally, we apply our model to the case of TiO in HD 209458b and show that the 
day-night cold trap would deplete TiO if it condenses into particles bigger than a few microns on the planet's night side, making it 
unable to create the observed stratosphere of the planet. 

Key words. Planets and satellites: atmospheres - Methods: numerical - Diffusion 



K>" 1. Introduction star (see Kirkpatrick (2005) for a review). As brown dwarfs pro- 
duce their own light, high-quality spectra of their atmosphere 
1988 was marked by the discovery of the first substellar ob- can be measured. Numerous atmospheric models have been built 
^ ject outside our solar system (Becklin & Zuckerman 1988). This ^o fit the data, from simple models with a reduced set of free pa- 
- - discoveiy was followed by numerous other discoveries of cool rameters (Ackerman & Marley 2001) to sophisticated models in- 
brown dwarfs, whose spectra differ significantly from stars. This eluding the detailed physics of condensation, growth and settling 
led to the definition of three new spectral classes beyond the stel- particles (Woitke & Helling 2004) 

lar classes O to M: the L T and Y dwarfs, ranging from temper- ^j^^ discovery and characterization of exoplanets followed 
atures of 2300 K to 350 K (Ctishing et al. 201 1). To obtain an ^j^^^ ^^j^^^^^ ^j^^ ^^^^ ^^^^^^^^^^ ^ ^^^^^ ^^^^^^ ^ ^^.^ 
understanding of these brown dwarfs and the mechanisms that ^^^^^ ^ ^ ^ q^^j^^ ^^995) opened the trail for the 
cause the transition from one spectral class to another, it is fun- discovery of hundreds of exoplanets. Some years later, atmo- 
damental to take into account condensation processes: m cold ^^^^.^ eharacterization of these objects became available, both 
atmospheres molecules can form, then condense and rain out, ^^^^ ^^^^-^ spectroscopy (Charbonneau et al. 2002) and from 
impacting the observed spectral features of these objects. In par- ^^^^^ detection of the planet's thermal emission (Charbonneau 
ticular, the M to L transition is marked by the disappearance of 2005). Although the global-mean effective temperatures of 
the titanium oxide TiO bands, which is understood by its trans- ^j^^^^ j^^^^^^ ^^^^^ ^^^^^ ^^^^^^ ^ 
formation into T1O2 and condensation into perovskite CaTi03 f^^ence is that hot Jupiters are strongly irradiated. Depending on 
(Lodders 2002). The L to T transition is understood as the switch .^^.^^^^ ^^^jj^^. ^^^^ ^^^emy rotation rate, and other factors, 
from a cloudy atmosphere to a cloud-free atmosphere as the atmospheric circulation models show that this day-night heat- 
cloud layer migrates mto the deep, unobservable regions of the -^^ ^^^^-^^^ j^^^^ pressures, to nightside temperatures 
that are at least ~1000 K colder than dayside temperatures (e.g. 



Send offprint requests to: V. Parmentier Showman et al. 2008, 2009; Dobbs-Dixon et al. 2010; Rauscher 



1 



Pamientier et al.: 3D mixing in hot Jupiter atmospheres I 



& Menou 2010, 2012c; Heng et al. 201 la,b; Pema et al. 2012). 
As a result, a wide variety of condensable species that are sta- 
ble in the gas phase on the dayside may condense on the night- 
side, leading to the formation of particles there. In the absence 
of atmospheric vertical mixing, such particles would gravitation- 
ally settle, depleting the atmosphere of these species on both the 
dayside and the nightside. Sufficiently strong vertical mixing, 
however, may keep these particles suspended in the atmosphere, 
allowing them to sublimate into gas phase in any air transported 
from nightside to dayside. Thus, the existence of this nightside 
"cold trap" is crucial for understanding not only the existence 
of hazes on hot Jupiters (e.g. Pont et al. 2012) but also the gas- 
phase composition of the atmosphere on both the dayside and 
the nightside. 

These arguments are relevant to a wide range of titanium 
and vanadium oxides, siUcate oxides, and other species. In par- 
ticular, chemical-equilibrium calculations show that, at temper- 
atures of ~ 1000-2000 K, there exist a wealth of condensates in- 
cluding Ti203, TisOs, and Ti407 (among other titanium oxides), 
MgAl204, Mg2Si04, MgSiOs, NaAlSisOg, KAlSisOg, and sev- 
eral phosphorus oxides (Burrows & Sharp 1999; Lodders 2002). 
Understanding the possible existence of gas-phase TiO, Na, K, 
and other species on the dayside would therefore appear to re- 
quire an understanding of the nightside cold trap. 

A particularly interesting problem in this regard is the 
existence of exoplanet stratospheres. For some transiting hot 
Jupiters, Spitzer IRAC secondary-eclipse observations indicate 
the presence of thermal inversions (stratospheres) on these plan- 
ets' daysides (Knutson et al. 2008). These stratospheres are gen- 
erally thought to result from absorption of starlight by strong 
visible/ultraviolet absorbers, but debate exists about the specific 
chemical species that are responsible. Hubeny et al. (2003) and 
Fortney et al. (2008) showed that, because of their enormous 
opacities at visible wavelengths, the presence of gaseous tita- 
nium and vanadium oxides can lead to stratospheres analogous 
to those inferred on hot Jupiters. 

However, a crucial question is whether the nightside cold 
trap would deplete the atmosphere of TiO, preventing this 
species from serving as the necessary absorber (Showman et al. 
2009; Spiegel et al. 2009). 

As pointed out by several authors, there exists another possi- 
ble cold trap. ID radiative-transfer models suggest that, on some 
hot Jupiters, the global-mean temperature-pressure profile be- 
comes sufficiently cold for condensation of gaseous TiO to oc- 
cur at pressures of tens to hundreds of bars (e.g., Fortney et al. 
2008). Even though a stratosphere on such a planet would be 
sufficiently hot for TiO — if present — to exist in the gas phase, 
the condensation of TiO and downward settling of the resulting 
grains at ~ 1 0- 1 00 bars might prevent the existence of TiO in the 
atmosphere (Showman et al. 2009; Spiegel et al. 2009). 

The strength of these two cold traps is given by a compe- 
tition between gravitational settling and upward mixing. The 
vertical mixing results from complex 3D flows. When using a 
ID-model, the vertical mixing is usually considered to be diffu- 
sive only, and parametrized by a diffusion coefficient K^^. Several 
studies give an estimate for this vertical diffusion coefficient in 
hot- Jupiter atmospheres. Heng et al. (201 lb) uses the magnitude 
of the Eulerian mean streamfunction as a proxy for the strength 
of the vertical motions and derived a vertical mixing coefficient 

of the order of » 10^ m Is. Some other works about chem- 
istry in hot- Jupiters decided to use an estimate based on the root 
mean squared vertical velocity either the local value (Cooper 
& Showman 2006) or the planet-averaged value (Moses et al. 



2011). However, these estimates are crude, and there is a need 
for theoretical work to more rigorously characterize the vertical 
mixing rates in hot Jupiter atmospheres. 

In this study we are interested in the horizontal cold trap. We 
we present three-dimensional general circulation model (GCM) 
experiments of HD 209458b to model a chemical species that 
condenses and settles on the nightside of the planet. We show 
that mixing in hot- Jupiter atmospheres is dominated by large- 
scale circulation flows resolved by the GCM. Finally, from the 
3D model, we derive the values of an effective vertical mixing 
coefficient, representing the global mixing of the planet. We also 
compare these 3D models to an idealized ID model parameter- 
izing mixing using an eddy diffusivity, with the goal of estimat- 
ing an effective eddy diffusivity for the mixing rates in the 3D 
models. These are the first circulation models of hot Jupiters to 
include a representation of clouds. 

2. 3D model 

Here, we use a state-of-the-art 3D circulation model, coupled to 
a passive tracer representing a condensable species, to determine 
how the interplay between dynamical mixing and vertical set- 
tling controls the spatial distribution of condensable species on 
hot Jupiters. Although the day/night cold trap should be present 
in most hot Jupiters, for concreteness, we must select a particu- 
lar system to investigate. HD 209458b is among the best studied 
hot Jupiters. Sitting at the boundary between pM and pL planets, 
it is believed to harbor a stratosphere (Knutson et al. 2008) and 
a strong day-night temperature contrast (Showman et al. 2009). 
We decided to use it as our reference model in this study, keeping 
in mind that most of the mechanisms discussed here should ap- 
ply to all hot Jupiters. To model the atmosphere of HD 209458b 
we use the 3D Substellar and Planetary Atmospheric Radiation 
and Circulation (SPARC/MITgcm) model of Showman et al. 
(2009), which couples the plane-parallel, multi-stream radia- 
tive transfer model of Marley & McKay (1999) to the MITgcm 
(Adcroft et al. 2004). 

2.1. Dynamics 

To model the dynamics of the planet we solve the global, three- 
dimensional primitive equations in spherical geometry using 
the MITgcm, a global circulation model for atmosphere and 
oceans developed and maintained at the Massachusetts Institute 
of Technology. The primitive equations are the standard equa- 
tions used in stably stratified flows where the horizontal dimen- 
sions greatiy exceed the vertical ones. In hot Jupiters, the hor- 
izontal scales are 10^ - lO^m whereas the vertical scale height 
fall between 200 and 500km leading to an aspect ration of 20 
to 500. In order to minimize the constraints on the timestep by 
the CFL criterion, we solve the equations on the cubed-sphere 
grid as described in Adcroft et al. (2004). The simulations do 
not contain any explicit viscosity nor diffusivity. However, in or- 
der to smooth the grid noise and ensure the stability of the code 
we use a horizontal fourth-order Shapiro filter (Shapiro 1970). In 
the vertical direction, no filtering process is applied. For a more 
detailed description of the equations see Showman et al. (2009). 

We use a gravity of 9.81 ms"^, a planetary radius of 
9.44 X 10^ m and a rotation rate of 2.06 x 10"^ s"' (implying 
a rotation period of 3.5 days). The average pressure ranges from 
200 bars at the bottom of the atmosphere to ptop at the top of the 
second-highest level, with the uppermost level extending from 
a pressure of /^top to zero. In most models, ptop is 2 fO-bar with 
53 vertical levels. In some models — ^particularly those with the 
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largest cloud particle size — we adopt ptop of 20i^bar with 47 ver- 
tical levels. In either case, this leads to a resolution of almost 
three levels per scale height. We use a horizontal resolution of 
C32, equivalent to an approximate resolution of 128 cells in lon- 
gitude and 64 in latitude and a timestep of 15 s. We reran some 
models at C64 resolution (equivalent to an approximate resolu- 
tion of 256x 128 in longitude and latitude) to check convergence. 

2.2. Radiative transfer 

The radiative transport of energy is calculated with the plane- 
parallel radiative transfer code of Marley & McKay (1999). The 
code was first developed for Titan's atmosphere (McKay et al. 
1989) and since then has been extensively used for the study of 
giant planets (Marley et al. 1996), brown dwarfs (Marley et al. 
2002; Burrows et al. 1997), and hot Jupiters (e.g. Fortney et al. 
2005, 2008; Showman et al. 2009). We use the opacities de- 
veloped by Freedman et al. (2008), including more recent up- 
dates, and the molecular abundances described by Lodders & 
Fegley (2002) and Visscher et al. (2006). As in the models of 
HD 209458b presented by Showman et al. (2009) and Showman 
et al. (2012), our opacity tables include gas-phase TiO and VO 
whenever temperatures are locally high enough for TiO to re- 
side in gas phase. Note that, because of the assumption of local 
chemical equilibrium, we do not consider the effect of cold traps 
on the atmospheric composition, and thus on the opacities. TiO 
opacities are always taken into account where temperatures are 
high enough for TiO to exist in the gas phase. This causes a warm 
stratosphere on the dayside of our modeled planet, regardless of 
the 3D distribution of our tracers (to be described below). 

We treat the opacities using the correlated-A; method (e.g.. 
Goody 1961). We consider 11 frequency bins for the opacities 
ranging from 0.26 to 300 |J.m; within each bin, opacity infor- 
mation from typically 10,000 to 100,000 frequency intervals is 
represented statistically. Note that this ^-coefficient approach is 
much more accurate than other multi-bin radiative transfer ap- 
proaches recently adopted in the literature (Dobbs-Dixon et al. 
2012; Dobbs-Dixon & Agol 2012). Showman et al. (2009) pro- 
vide a detailed description of the radiative-transfer model and its 
implementation in the GCM. 

2.3. Tracer fields 

Our target species is represented by a passive tracer field. Tracer 
fields are often used in GCMs to follow the concentration of 
a chemical species such as water vapor or cloud amount in the 
atmosphere or salinity in the oceans. Cooper & Showman (2006) 
were the first to include a chemical tracer in a circulation model 
of hot Jupiters, in their case to investigate the quenching of CO 
and CH4 due to atmospheric mixing. Our study represents the 
first hot- Jupiter circulation model to include a representation of 
clouds. 

The tracer field is advected by the flow calculated in the 
GCM. Thus the tracer abundance x is given by the continuity 
equation: 

Dy 

-^=S (1) 
Dt 

Here, ^ represents the mole fraction of the tracer, i.e., the number 
of molecules of the species in a given volume, either in gaseous 
phase or trapped in condensed particles, with respect to the total 
number of atmospheric molecules in that volume. For simplic- 
ity the value of x is normalized to its initial value in the deep 
layers of the planet. In the equation, 5 is a source term and the 



total derivative is defined by D/Dt = d/dt + v • V + where 
t is time, v is the horizontal velocity, V is the horizontal gradi- 
ent operator on the sphere, w = Dp/Dt is the vertical velocity 
in pressure coordinates, and p is pressure. We model a condens- 
able species which is gaseous on the dayside of the planet, and 
trapped into particles of size a on the nightside. Thus the source 
term S represents the gravitational settling of these particles and 
is given by 



_ 



on the dayside 
-g-^ — on the nightside 



(2) 



where z is the height, increasing upward, p is the density of the 
air and Vf is the settling velocity of the particles defined by 
equation 3. This velocity depends on the size of the particles, 
which is determined by the complex microphysics of conden- 
sation, out of the scope of this study. Thus we treat a as a free 
parameter in our model. We model spherical particles of radii 
0.1, 0.5, 1, 2.5, 5 and 10 \im. Equation 2 describes a simple, 
bimodal mechanism for the condensation of chemical species 
on the nightside. Although this scheme is highly simplified (ig- 
noring the detailed temperature- and pressure-dependence of the 
condensation curves of possible condensates in hot Jupiter at- 
mospheres), it reflects the fact that a wide range of species will 
reside in gaseous form on the dayside yet condensed form at low 
pressures on the nightside. 

The tracers are considered passive, which means they do not 
influence either the dynamics or the radiative transfer of the sim- 
ulation. This ignores a role for possible radiative feedback mech- 
anisms - discussed in Sect. 6 - but represents a necessary first 
step toward understanding how dynamics controls the 3D distri- 
bution of a condensable species. 

2.4. Settling velocity 

We assume that the target chemical species condenses into 
spherical particles of radius a that reach inmiediately their ter- 
minal fall speed, which is given by (Pruppacher & Klett 1978): 



Vf 



2j8a^g(Pp-p) 
9t] 



(3) 



where is the viscosity of the gas, g is the gravitational accel- 
eration of the planet, pp is the density of the particle and p the 
density of the atmosphere. Vf is positive when the particles goes 
downward. The Cunningham slip factor, jS, accounts for gas ki- 
netic eff'ects that become relevant when the mean free path of the 
atmospheric molecules is bigger than the size of the falling parti- 
cle. This factor has been measured experimentally by numerous 
experiments. We adopt the expression from Li & Wang (2003) 
as done by Spiegel et al. (2009). 



P=l+ Kt^(l.256 + 0.4e-''/^'^) 



(4) 



where the Knudsen number ^Tn is the ratio of the mean free path 
to the size of the particle : 



(5) 



For a perfect gas, the mean free path can be expressed as 
(Chapman & Cowling 1970): 



ksT 1 



(6) 
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a (jtm) 



Fig. 1: Terminal velocity of a falling particle as a function of 
pressure and particle size for an H2-atmosphere at 1000 Kand a 
value of the gravitational acceleration of 15 m s"^. 

with d the diameter of the gas molecules, P and T the pressure 
and temperature of the gas and kg the Boltzmann constant. 

In the limit of a high-density atmosphere, yS — > 1 and the 
terminal speed Vj becomes the Stokes velocity Vstokes- 

For low density gases, the dynamical viscosity is indepen- 
dent of pressure and can be expressed as a power law of the 
local temperature with an exponent varying between 1/2 in the 
hard-sphere model to near unity, depending on the strength of 
the interactions between the molecules. Following Ackerman & 
Marley (2001), we use the analytical formula given by Rosner 
(2000) for the viscosity of hydrogen : 

5 yl^^^iksTlef '^ 

with d the molecular diameter, m the molecular mass and e 
the depth of the Lennard-Jones potential well (for H2 we use 
2.827 X lO"'" m and 59.7kB K respectively). The power law be- 
havior of the viscosity remains valid for temperatures ranging 
from 300 K to 3000 K and for pressures less than 100 bars (Stiel 
& Thodos 1963). At higher temperature, ionization of hydrogen 
becomes relevant and the viscosity reaches a plateau. However 
the temperatures of the model are everywhere less than 3000 K 
and so equation 7 remains valid. 

Figure 1 displays the resulting terminal velocity as a function 
of pressure and particle size. Two different regimes are observed. 
For Knudsen numbers smaller than unity, the terminal velocity 
is independent of pressure, whereas for Knudsen numbers ex- 
ceeding unity, the terminal velocity is inversely proportional to 
the pressure. At low pressure and for particles bigger than a few 
tens of micrometers, the Reynolds number becomes higher than 
unity, and Equation 3 is no longer valid. However, we show in 
Appendix A that these differences remain smaller than one order 
of magnitude and confined to a small parameter space thus we 
decided to neglect them for this work. 

2.5. Integration time — limitation of the study 

A challenge for any 3D numerical integrations of hot- Jupiter dy- 
namics is the wide range of timescales exhibited by these at- 
mospheres. This is true for the radiative time constant, which 




0.1 0.5 1 2.5 5 10 100 



a (|jm) 

Fig. 2: Settling timescale as a function of pressure and particle 
size for the same conditions as in figure 1 . Also shown are the 
advective timescale (black curve) and the simulation timescale 
(blue curve). In the shaded region (below Ibar), the tracers are 
considered gaseous and the settling timescale is not relevant for 
our study. 



varies significantly from low to high pressure (Iro et al. 2005; 
Showman et al. 2008). Nevertheless, Showman et al. (2009) 
showed that the computed light-curve of HD 189733b changes 
little for integration time longer than hundreds of days, indicat- 
ing that the dynamics in the millibar regime has stabilized. As 
we integrate ~ 1400 days, we consider the dynamics of the planet 
to be spun up at pressures less than ~100 mbar. For the settling 
of particles, another timescale must be taken into account. The 
particle settling timescale can be defined as the time for the par- 
ticles to fall one atmospheric scale height. To obtain a correct 
picture of the problem at a given location of the planet, we need 
to integrate the simulation for at least several times longer than 
the settling timescale at this location. As we can see in fig. 2 the 
settling timescale ranges from tens of seconds for big particles at 
low pressure to tens of years for small particles at high pressure. 
Due to computational limitation it is not possible to run the sim- 
ulation long enough to ensure that every considered tracer field 
has reached a statistical steady state for the full range of particle 
sizes we consider. The integration during 1400 days allows to 
calculate the steady state for every particle size at pressure lower 
than ten millibar and at every pressure for particle size bigger 
than 2 \im. 

We define the advective timescale as the time for a parcel 
in the main jet stream to cross one hemisphere of the planet. If 
the advective timescale exceeds the settling timescale at a given 
level, the particles at these levels will fall several scale heights 
while on the nightside. We thus expect that these levels will 
be depleted. Conversely when the advective timescale is shorter 
than the settling timescale, the coupling between the flow and the 
particle is essential and particles' behavior cannot be predicted 
easily. Typical advective times in our models are ~24 hours; this 
is marked by the upper thick black curve in fig. 2. We expect 
depletion to occur above this line. 
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Fig. 3: Temperature (left panel, colorscale), vertical velocities (right panel, colorscale) and horizontal winds (arrows) in our model 
of HD 209458b. Positive velocities are upward. The top three panels show the flow at three different pressure (O.lmbar, Imbar and 
lOmbar). The bottom panel show the vertical velocities versus longitude and pressure along the equator. The substellar point is at 
longitude, latitude (0°, 0°), the dayside is between -90° and +90°. All the quantities are time averaged. 



3. Results 

3.1. Dynamical regime 

We run the simulations for 1400 days and calculate the time 
average for all the variables over the last 400 days only, once 
the simulation reaches a statistical steady state at upper levels. 
Understanding the flow structure is essential to understanding 



the Lagrangian advection of particles, thus we first present a 
brief description of the dynamics. A more complete description 
of the circulation in our simulations is presented in Showman 
et al. (2009). The temperature structure is shown in fig. 3. A 
hot stratosphere is visible at altitudes above the 10-mbar level, 
due to the strong visible-wavelength absorption by titanium ox- 
ide present in solar abundances in the simulation. Temperatures 
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reach ~2200 K at low pressures near the substellar point. By 
contrast, the temperatures deeper than 10 mbar in the dayside 
are relatively temperate. This could lead to the presence of a ver- 
tical cold trap, not considered in this study. The day-night tem- 
perature contrast becomes significant at pressures less than ~100 
mbar, reaching 1600 K near the top of the model. This large tem- 
perature difference results from the short radiative timescale at 
low pressures in comparison to dynamical timescales (Iro et al. 
2005; Cooper & Showman 2005; Showman et al. 2008). 

The horizontal flow on isobars comprises an eastward (super- 
rotating) jet close to the equator and a day-to-night flow pattern 
at higher latitude. As the pressure increases, the radiative time 
constant increases and the jet extends to higher latitudes. The 
mean vertical velocities exhibit planet-wide variations (fig. 3). 
The highest velocities coincide with strong horizontal conver- 
gence of the flow and occur mostly at the equator West of the 
anti-stellar point, the convergence of the day-to-night circulation 
forces strong downwelling motions. As described in Rauscher 
& Menou (2010), this convergence point appears in a range of 
hot Jupiter GCMs of varying complexity. It is usually associated 
with a shock-like feature (Heng 2012). Our simulations assume 
local hydrostatic equilibrium and thus can not treat properly the 
physics of shocks. To date, no model of hot Jupiter atmospheric 
dynamics can handle shocks properly. Yet, a similar wind con- 
vergence pattern appears when considering the non-hydrostatic 
nature of the flow as can be seen in Dobbs-Dixon et al. (2010). 
Between the substellar point and the west terminator, there are 
additional convergence/divergence points associated with the jet, 
leading to a region of strong ascending motion ~40° of longitude 
west of the substellar point, and a broad region of descending 
motion west of that. These vertical flows remain coherent over 
several orders of magnitude in pressure, giving them the poten- 
tial to transport vertically large quantities of material. Outside 
of these points of strong vertical motions, the vertical velocities 
are more than one order of magnitude smaller They are mostly 
upward on the dayside and downward on the nightside. 

3.2. Spatial distribution of condensabie species 

Our simulations show that, as expected, condensation and parti- 
cle settling on the nightside depletes the tracer from upper levels 
relative to the abundances at depth — an eff'ect that is stronger 
for larger particles. This is illustrated in fig. 4 (solid curves), 
which shows the global-mean tracer abundance (averaged hori- 
zontally on isobars) versus pressure for simulations with particle 
sizes ranging from 0.1 fim to 10/im. In all cases, the horizontally 
averaged tracer abundance decreases with altitude. The deple- 
tion is modest for the smallest particle size (0.1 fim), but for the 
largest particle sizes, tracer abundances at the top are almost two 
orders of magnitude smaller than abundances at the bottom. A 
useful metric is the "50% depletion pressure," that is, the pres- 
sure above which the global-mean tracer abundance is less than 
50% of the deep abundance. This pressure is only 10/ibar for a 
particle size of 0.5 //m but is 0. 1 bar for a particle size of 10 //m. 
We also note that, for all the models shown in fig. 4, the particles 
settling times near the top of the model are much less than our 
integration times; at low pressures, the tracer abundances have 
reached a statistical equilibrium where downward transport of 
tracer due to particle settling is balanced by upward mixing of 
tracer by the large-scale dynamics. 

The tracer abundance on isobars exhibits a strong spatial 
variation as can be seen in fig. 5. Although a day/night pattern 
is imposed in the tracer source/sink (with particle settling on the 
nightside but not the dayside), the three-dimensional advection 




Mean tracer abundance 



Fig. 4: Planet-wide, time averaged tracer abundances on isobars. 
The abundances are normalized to the abundance in the deepest 
layer We compare the value of the 3D model (solid lines) and 
the fit using the ID model (dashed lines). 



of the tracer by the atmospheric winds leads to a complex tracer 
distribution that does not exhibit an obvious day-night geometry. 
The main pattern appears to be an equator-to-pole gradient, with 
large zonal-mean abundances at the poles, and smaller zonal- 
mean abundances at the equator This is particularly true around 
O.lmbar. Significant longitudinal tracer variations also occur; at 
~1 mbar (fig. 5), these variations are particularly prominent at 
high latitudes. Interestingly, these variations are phase shifted in 
longitude relative to the day-night pattern with maximum (min- 
imum) peak tracer abundances occurring ~60-80° of longitude 
east of the substellar (antistellar) point. 

On top of these main patterns, we clearly see two points de- 
pleted in tracers at the equator These two points correspond to 
points of horizontal convergence of the flow and high down- 
welling motions as discussed in the previous section. This corre- 
lation between strong downwelling motions and low tracer abun- 
dances arises naturally in the presence of a background vertical 
gradient of tracer abundances. Due to their settling on the night- 
side, the local tracer abundance generally decreases with height 
and thus any downwelling motion would carry a parcels of gas 
depleted in tracer whereas any upwelling motion should carry 
parcel of gas enhanced in tracer 

3.3. Geometry of the mixing 

As discussed in the previous section, the tracer abundance is 
not null everywhere in the planet, which implies that tracers did 
not rain out during the simulation time. Yet the integration time 
greatly exceeds the fall times at low pressure for all particle sizes 
considered, and everywhere throughout the domain for particle 
sizes exceeding a few yum. Therefore vertical mixing must hap- 
pen in order to keep these particles lofted. This vertical mixing 
is characterized by an upward dynamical flux of tracers that b al- 
ances the downward flux due to the gravitational settling in the 
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Fig. 5: Time averaged tracer abundance at two different pressures and for different particle sizes. 
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nightside. The upward dynamical flux of tracers across isobars 
can be calculated as (-^) where u) is the vertical velocity in 
pressure coordinates, x the tracer abundance, g the gravity of the 
planet and the brackets denote the horizontal (global) average on 
isobars. As can be seen in fig. 6, the upward flux of tracers due 
to the dynamics (solid lines) balances nicely the downward flux 
due to settling (dashed lines), showing that the simulation did 
reached a quasi steady-state. 



0.1^m : 

0.5M.m i 
1|j.m : 

2.5^im J 

5|j,m 




-3-2-10123 
Fluxes accros isobars (kg/s) - negative is upward ^rfi 



Fig. 6: Time averaged fluxes of tracers across isobars. Upwelling 
fluxes due to the dynamics (solid lines) balance the downwelling 
fluxes due to the gravitational settling of the tracers (dashed 
lines). Negative values are upward fluxes. 

Hot Jupiter atmospheres are heated from above and thus are 
believed to be stably stratified. Then, vertical mixing can not 
be driven by small scale convection as it is the case in the deep 
atmosphere of brown dwarfs (Freytag et al. 2010) and the gas gi- 
ants of the solar system. As explained in Sect. 2. 1, we don't have 
any parametrization of sub-grid scale mixing. Thus we do not ac- 
count for mixing induced by small-scale turbulence and gravity 
wave breaking, two mechanisms that are believed to dominate 
the mixing in the radiative part of brown dwarfs atmospheres 
(Freytag et al. 2010). Rather, the upward flux of tracer in our 
model is due to the large-scale, resolved flow of the simulation. 

Given that mass is conserved, any upward flux of gas is com- 
pensated by a downward flux of gas. Thus, if the tracer concen- 
tration were horizontally homogeneous on isobars, there would 
be no net upward flux of tracer through that isobar For a net 
upward flux of tracers across isobars to occur, there must be a 
correlation between the horizontal distribution of the tracers and 
the vertical velocities. Dynamics will produce an upward flux 
if — on an isobar — ascending regions exhibit greater tracer abun- 
dance than descending regions. In other words, an upward tracer 
flux due to dynamics will occur only if 0^ - {)()) to < Q where 
10 is the upward velocity in pressure coordinates, x is the tracer 
abundance and the brackets are the mean over one isobar (note 
that negative oj implies upward motion). Given a vertical gra- 
dient of (x) such that the abundance of {x) decreases upward, 
an upward flux of gas will naturally bring enhanced material 
whereas a downward flux will naturally advect parcels of gas 
depleted in tracers, thereby creating the correlation between co 
and X favorable for upward tracer transport. Figure 7 shows the 



relative contribution to the upward mixing versus longitude and 
latitude at a given isobar, which can be estimated by the quantity 

P^'^iXzMl (8) 

Mass conservation in the primitive equations equation implies 
that {(jj{x)) - 0: the advection of (^) does not contribute to the 
net (horizontally averaged) upward flux of material and so we re- 
move the contribution of this term when defining equation 8. The 
quantity is the mean upward flux of material across iso- 

bars, thus the quantity F represents the local contribution to the 
global upward flux on isobars. It is normalized such that {F) - 1 . 

The strength of the mixing varies significantly with longitude 
and latitude. Both upward and downward fluxes are one order 
of magnitude larger than typical values in a handful of specific 
small areas across the planet — particularly at the two points of 
horizontal convergence and strong vertical velocities described 
in Sect. 3.1. This vertical flows remains coherent over several 
order of magnitude in pressure (see fig. 3), acting like a vertical 
"chimney" where eflicient transport of material can be achieved. 

In summary, the mechanism by which the large-scale, re- 
solved atmospheric circulation transports tracer upward is ex- 
tremely simple and straightforward. The settling of particles 
leads to a mean vertical gradient of tracer abundance, with, on 
average, small tracer mixing ratios aloft and large tracer mix- 
ing ratios at depth. Given this background gradient, advection 
by vertical atmospheric motions — whatever their geometry — 
autmomatically produces a correlation between co and (x - {x)) 
on isobars, with ascending regions exhibiting larger values of 
ix - {x)) than descending regions. In turn, this correlation auto- 
matically causes an upward dynamical net flux of tracer when 
averaged globally over isobars. In statistical steady state, this 
upward dynamical flux balances the downward transport due to 
particle settling and allows the atmospheric tracer abundance to 
equilibrate at finite (non-zero) values despite the eff'ect of par- 
ticle settling. The mechanism does not require convection, and 
indeed, the vertical motions that cause the upward transport in 
our models are resolved, large-scale motions in the stably strat- 
ified atmosphere. These vertical motions are a key aspect of 
the global-scale atmospheric circulation driven by the day-night 
heating contrast. 

3.4. Time variability 

Besides the spatial variability at a given time, the model exhibit 
significant temporal variability, both in the 3D flow and espe- 
cially in the tracer field. The equatorial jet exhibits a global os- 
cillation pattern as can be seen in fig. 9. At the convergence point 
west of the substellar point, the jet orientation can be toward the 
north, the south or well centered on the equator Whereas Dobbs- 
Dixon et al. (2010) described a variation in longitude with time 
of the convergence point of the flow in the nightside, we see a 
variation in latitude of this convergence point and interpret it as 
a result of the larger oscillation of the jet itself. 

The tracer abundances at specific locations on the planet 
exhibit strong temporal variability. This is illustrated in fig. 9, 
which shows the tracer abundance at 0.1 and 1 mbar over the 
globe at several snapshots in time for a model where the ra- 
dius of particles on the nightside is 2.5 fim. Significant varia- 
tions in tracer abundance are advected by the equatorial jet and, 
at high latitudes, by the day-to-night flow, leading to large lo- 
cal variations in time. In many cases, the strongest tracer vari- 
ability seems to involve regional-scale structures with typical 
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Fig. 7: 2D distribution of the mixing efficiency on the Imbar iso- 
bar. We plot F - ""^^^"1^^ ■ Positive is either an upward flux of 
gas enhanced in tracers compared to the horizontal mean or a 
downward flux of gas depleted in tracer compared to the hor- 
izontal mean value. Negative is either a downward flux of gas 
enhanced in tracer or an upward flux of gas depleted in tracer 
Thus a positive value enhances the (horizontally averaged) up- 
ward flux, whereas a negative value diminishes it. Theses fluxes 
are normalized to the mean upward flux of tracer and the global 
mean on an isobar of the plotted quantity is 1 . 



sizes of ~l-3xlO'*km but also includes hemispheric-scale fluc- 
tuations (e.g., in the abundance averaged over the day or night) 
and between the northern and southern hemispheres. Around the 
substellar point, the tracer abundance can vary by up to 50%, 
whereas along the terminator, this temporal variation can reach 
75% relative to the mean value. Such variability — if it occurs in 
radiatively active species like TiO — has important implications 
for secondary-eclipse and transit observations, which probe the 
dayside and terminator, respectively. 

Figure 8 sheds light on the different timescales at which this 
variability occurs. The top panel shows the tracer abundance av- 
eraged vertically between 0. 1 and 1 mbar and horizontally over 
a circular patch with a radius of 45° centered on the substellar 
point; this gives a sense of how the tracer abundance would vary 
in secondary-eclipse measurements probing the dayside. The 
bottom panel shows the tracer abundance averaged vertically be- 
tween 0.1-1 mbar and horizontally around the terminator, in- 
cluding all regions within 5° of the terminator itself; this gives a 
sense of how the tracer might vary in transit measurements. The 
variability exhibits two characteristic timescales: a short (fast) 
timescale of order of days and a long (slower) timescale of ~50 
to 100 days. The bigger the particles on the nightside, the bigger 
the amplitude of the variations. In these models, the amplitude of 
the long-period variations exceeds those of the short-period vari- 
ations by a factor of ~2-3. The long-timescale variations exhibit 
similar amplitudes in the dayside and terminator time series. The 
short-period oscillations exhibit stronger amplitude at the termi- 
nator and seem related to variations of the flow itself, such as 
the oscillation of the jet described previously. Figure 8 suggests 
that radiatively active tracer species that can condense on the 
nightside, such as TiO or silicates, could lead to detectable time 
variations in transit or secondary eclipse spectra. The amplitude 
of this variability will depend on the type of tracer being consid- 
ered (see Sect. 6) and may vary from planet to planet depending 
on the availability of the considered species. However we can 
predict the expected period of these variations; some days for 
the small amplitude ones and fifty to one hundred days for the 
biggest ones. 
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Fig. 8: Time dependence of the mean tracer abundance averaged 
between Imbar and 0.1 mbar and averaged horizontally over a 
circular region of radius 45° centered on the substellar point 
(top) and averaged horizontally around the terminator, includ- 
ing all regions within ±5° of the terminator (bottom). The top 
panel is relevant for inferring the presence of a stratosphere (see 
Sect. 6.1), albedo variations (see Sect. 6.2) or secondary eclipse 
measurements (see Sect. 6.3). The bottom panel is relevant for 
transit spectroscopy measurements. 



3.5. Limb profile 

Transit observations are sensitive to atmospheric composition 
near the terminator, and it is therefore of interest to characterize 
the distribution of our tracer species at the terminator The possi- 
bility of variations in chemical composition between the leading 
and trailing limbs (as seen during transit) has been discussed in 
a variety of studies (e.g., Iro et al. 2005; Fortney et al. 2010). 
However, these studies did not investigate the particular deple- 
tion of species due to the interaction between their condensation 
and the atmospheric dynamics. Our model allows the first quan- 
titative estimate of how dynamics and nightside condensation 
affect the abundance of condensable species at the terminator. 

Figure 1 1 shows the tracer abundances at the terminator at a 
snapshot in time for our models with particle sizes of 0.1, 0.5, 1, 
2.5, 5, and IQ^m. Angle represents angle around the terminator 
and the radial coordinate represents log pressure. In agreement 
with fig. 4, the tracers tend to be depleted from upper levels, par- 
ticularly in models where the particles on the nightside are larger 
Moreover, fig. 1 1 demonstrates that significant spatial variations 
occur along the terminator Depletion occurs first at the equator 
along the leading limb, corresponding to the terminator 90° west 
of the substellar point. The superrotating jet canies air depleted 
in tracer from the nightside directly to this region of the termina- 
tor, explaining why abundances are particularly depleted there. 
In contrast, air along most of the remainder of the terminator has 
arrived from the dayside, where no particle settling occurs, so 
depletion is less strong — particularly for particle sizes < 1 yum. 
Once the particle size becomes sufficiently large, however, de- 
pletion occurs everywhere along the terminator at upper levels 
regardless of whether the air arrived there from the dayside or 
nightside. 

Considering now the depth dependence of the terminator 
abundances, our results suggest two different zones (see figs. 10 
and 11): 
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Fig. 9: Tracer abundance (colorscale) and winds (arrows) at two different pressure levels and 3 different times of the simulation for 
the case a = 2.5 [a.m. 



- At altitudes above the 1-mbar level, the tracer abundance is 
homogeneous over most of the limb except for the east (trail- 
ing) equatorial limb that is strongly depleted. 

- At altitudes below the 1-mbar level, the east and west equa- 
torial limb are rather homogeneous; however, the east/west 
dichotomy shifts to higher latitudes and the west limb above 
45° is more depleted than the equivalent region of the east 
limb. 

Moreover, it should be noted that due to the shift of the hot 
spot, the temperatures at the east limb are higher than at the west 
limb, thus a given species is more likely to be gaseous and de- 
tectable on the trailing limb than on the leading limb. 

4. 1 D model of the day-night cold trap 

Although hot Jupiter atmospheres are inherently three dimen- 
sional, ID models continue to play a useful role for understand- 
ing the vertical thermal and chemical structure of these atmo- 
spheres. In particular, many groups have explored the chem- 
istry of hot Jupiters using ID models in which the vertical 
mixing caused by the large-scale dynamics is parameterized by 
a specified eddy diffusivity (e.g. Spiegel et al. 2009; Zahnle 



et al. 2009a,b; Youdin & Mitchell 2010; Line et al. 2010, 2011; 
Madhusudhan & Seager 201 1). In these studies, the chosen eddy 
diffusivity is ad hoc, with no convincing theoretical support. 
Although the vertical mixing in our 3D models is not diffusive 
in any rigorous sense, there is merit in comparing the results of 
our 3D models with ID models parameterized by eddy diffu- 
sivity. This will allow us to make approximate estimates of the 
magnitudes of eddy diffusivity — in the context of a ID model — 
that produce similar horizontal-mean behavior as our 3D mod- 
els. Such estimates of eddy diffusivity should help to guide pa- 
rameter choices in ID chemical models like those cited above. 
A comparison between our 3D models and ID diffusive models 
will also allow us to investigate how the horizontal-mean tracer 
depletion relates to the amplitudes of spatial tracer variation. 

Therefore, in this section, we present a simple ID model, 
including particle settling, with atmospheric mixing represented 
as an eddy diffusivity. 

4. 1 . System studied 

The presence of a superrotating, eastward equatorial jet is a dom- 
inant dynamical feature of many 3D circulation models of hot 



10 



Parmentier et al.: 3D mixing in hot Jupiter atmospheres I 




Fig. 10: Limb profiles of our tracer field from 200bar to 1 p.bar as seen during transit. Black circles are situated at lOObar, lObar, 
Ibar, O.lbar, O.Olmbar, Imbar, O.lmbar, O.Olmbar and l|xbar. The north pole is on top and the leading limb on the right. The 
abundances are normalized to the abundance of the deepest level. 



11 



Parmentier et al.: 3D mixing in hot Jupiter atmospheres I 



10 



10 



10^ 



^ 10"^ 

(U -2 

§ 10 

in 

^ -1 

Q- 10 



10" 



10^ 



10^ 




0.2 



0.4 0.6 0.8 
Mean tracer abundance 



Fig. 1 1 : Tracer abundance along the day-to-night terminator 
(equator in black, poles in red) and along the night-to-day ter- 
minator (equator in green, poles in blue) for the case a - 2.5i_im. 



Jupiters. This superrotating jet was first predicted by Showman 
& Guillot (2002), and later emerges from almost all 3D simu- 
lations of hot Jupiter atmospheres (Cooper & Showman 2005; 
Showman & Cooper 2006; Showman et al. 2008, 2009, 2012; 
Dobbs-Dixon & Lin 2008,?; Menou & Rauscher 2009; Rauscher 
& Menou 2010, 2012b,a; Pema et al. 2010, 2012; Heng et al. 
2011a,b; Lewis et al. 2010; Kataria et al. 2012) including ours 
(see Sect. 3.1) and has been theoretically understood (Showman 
& Polvani 201 1). A shift of the hottest point of the planet east- 
ward from the substellar point has been directly observed in 
several exoplanets (Knutson et al. 2007, 2009, 2012; Crossfield 
et al. 2010) and interpreted as a direct consequence of this jet. 
Thus we believe that any study of the day/night cold trap in hot 
Jupiter atmospheres must account for this feature. 

To include the presence of this jet in our model, we choose 
as a study system a vertical column of gas homogeneously ad- 
vected around the equator by the superrotating jet. Such a col- 
umn is transported from day to night and from night to day with 
a period T^dv = — ^ where Tajy is the advective timescale, /?„ 
is the planetary radius and Mjet the equatorial jet velocity, r^dv is 
around 48h for HD 209458b. 

We focus on a hypothetical chemical species which is 
gaseous on the dayside and condenses on the nightside. This 
species freely diffuses with a vertical diffusion coefficient K^^ 
on both the dayside and the nightside. Because we envision 
the species as condensed on the nightside, we additionally in- 
clude downward settling via Stokes-Cunningham drift on the 
nightside — but not on the dayside. The model includes no hori- 
zontal dimensions. Rather, we envision that the model captures 
the vertical diffusion in a given column as it as advects from 
day to night, and we therefore include the day-night difference 
in settling (with settling on the nightside but not dayside) as a 
time dependence. We assume this chemical species to be a mi- 
nor constituent of a Ha-atmosphere. Thus we neglect the latent 
heat released during the condensation. 



4.2. 1 D diffusion equation 

As before, x the local mole fraction of the target chemical 
species, giving moles of tracer species (whether in gaseous or 
condensed form) to the total moles of air in a given volume. On 
the dayside, the molecules of the target chemical species can 
freely diffuse with a diffusion coefficient K^z, according to the 
equation 



with p the density of the atmosphere and z the vertical coordi- 
nate. 

On the nightside, the molecules of the target species are 
trapped into particles that both diffuse and settle with the ve- 
locity described in 2.4. Thus x follows the same equation as 9, 
plus a source term describing the settling: 



dx 1 d I^^ Jx\ ^ 1 d(fixVf) 
dt pdz\ dzj p dz 



(10) 



with Vf the settling velocity. 

Assuming hydrostatic balance, we can use the pressure P as 
the vertical coordinate and the equation becomes: 



dt 



1 djpxVf) 
g dP 



Using the perfect gas law for the atmosphere we get: 



dx__ d_ 
dt dP 



mg_ 
knT 



d(^PxVf) 
dP 



(11) 



(12) 



We can define the diffusive time scale as Td = and a reference 
free fall time scale Ts - H/Vs with H - ^ the atmospheric 
scale height and Vs the Stokes velocity (see 2.4). We note that Td 
is a reference time scale and is not equal to the effective free fall 
time scale for high Knudsen numbers. 
Equations 9 and 10 become: 



dt p dz\ " dz 



dt pdz\ dz 







1 dpPx 
Ts dP 



on the day side 



on the night side 



(13) 



4.3. Time-dependent solution 

In order to solve the system of equations 13 for xiP^O, we 
need two boundary conditions. We assume the species to be 
well mixed below f i = 1 bar with a molecular abundance 
XQ. At the top of the atmosphere, we assume that no molecule 
crosses the upper boundary (i.e.^\p=p^^^^ -0). Then we can 
solve the system with an implicit time stepping code using 
GNU octave', an open-source, free software equivalent to 
MATLAB. Assuming that the column of gas spends Tadv/2 = 
24 h in each hemisphere, we can reach a periodic behavior where 
the initial condition is forgotten. While on the dayside, there is 
no settling and, at upper levels, the tracer diffuses upward. Thus 
at a given pressure the molecular abundance increases with time. 
This is shown by the red curves in fig. 12. While on the night 
side, the particles both diffuse vertically and settle downward 



http://www.octave.org 



12 



Parmentier et al.: 3D mixing in hot Jupiter atmospheres I 



at their terminal velocity and thus, at at given pressure level, the 
abundance decreases with time. This is shown by the blue curves 
in fig. 12. In the upper atmosphere, the settling timescale and the 
diffusion timescale become smaller than the advective timescale. 
The particles have time to settle several scale heights during the 
time required for the air column to cross the dayside. Thus, at 
upper levels, the molecular abundance vary strongly throughout 
the diurnal cycle, and the vertical profiles vary widely around the 
mean value. 

4.4. Steady-State solutions 

Our system being forced periodically, there is no steady-state so- 
lution stricto sensu. However, the mean of the tracer abundance 
over one period should remain constant. We can thus integrate 
the system 13 over one period and get an equation for ^, the 
mean molecular abundance over one period. Then, once the pe- 
riodic state is reached, the mean over one period of the tracer 
abundance in our column of gas, is the same as the mean over 
longitude at a given time, (x) and we obtain : 



OP 2 Ts 



(14) 



where the factor ^ appears because the source term is only inte- 
grated during the night. To derive Equation 14 we assumed the 
mean tracer abundance on the nightside is close to the global- 
mean tracer abundance. This is true given than fir^ and Td are 
much bigger than Tadv This approximation breaks down at low 
pressure and our analytical solution diverges from the real one 
as can be seen by comparing the black and the green curves in 
fig. 12. However, the discrepency is small such that the analytic 
solution is still a good representation of the time-mean of the full 
numerical solution of the ID model everywhere deeper than the 
~ 10-jubar level. C is a constant coming from the integration over 
pressure. When P goes to 0, x goes to 0. Moreover we assume 
that does not go to infinity. Then the constant must be zero 



and we obtain ; 



dP 2 P Ts ^' 



(15) 



To simpUfy the problem, we neglect the transitional regime 
for p from equation 4 and use the expression : 



j8 = 1 + l.656Kt, 



(16) 



Then choosing a functional form K^^ - K-i^2.o(Pq/P)", where a is 
a constant, we can solve equation 15. For a and a 1, we 
get : 



. 1 TA^P" -P" 



exp 



2a T, p;; 

1 Ti,l.656k^T P"-'^ -P"^-^ 
2(ff-l)77 yl2„acf PI 



(17) 



For the particular case of a constant K^z (fl - 0) the formula 
becomes : 

And in the case where is inversely proportional to P (a = 1): 

1 '"do l.656k^T 

WheTexr is the abundance of (x) at P - P^ and Td,, = H/Kzzf, 
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Fig. 12: Tracer abundance in the advected column of gas from 
the ID model as a function of time. From bottom to top we used 

K,, = (10^5 X 10^ 10-**) VI bar/Pm'/s and a particle radius 
on the night side of 1 |^m. We plot one profile every hour The 
blue ones are on the night side, the red ones on the dayside. The 
black curve is the mean over one period. The green line is the 
analytical model (equation 17 using a - 1/2) and the dotted line 
is the solution for particles that would be constantly falling (both 
in the dayside and in the night side). 
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4.5. Diffusivities needed to keep tracer suspended 

Because particle settling acts to transport condensates down- 
ward, the tracer only exhibits significant abundances in the up- 
per regions of the atmosphere if the eddy diffusion coefficient 
exceeds some critical value, which depends on the particle size 
of the condensates. Here we solve for an approximate analytical 
expression for this critical magnitude as a function of particle 
size and other parameters. 

Assuming a constant vertical diffusion coefficient, we can 
use equation 18 to derive an expression for the K^z needed to 
achieve a given molecular abundance ;^'f at a given pressure P{. 
We first notice that equation 18 is composed of two terms. The 
first exponential is given by the Stokes regime (/? » 1) whereas 
the second term is given by the Cuningham regime (J3 » 1). As 
can be seen in fig. 1, particles smaller than 10 |j.m should be in 
the Cunhingham regime at f = 1 mbar. Thus the second term is 
dominant and we can use as a condition : 



1 Td l.656kBT 
2ts yJlnacP- 



P Pi) 



■ InOfiim). 



(20) 



Assuming Pi » P\im, replacing Td and Ts by their expres- 
sions, and using equation 7 for the viscosity and equation 6 for 
the mean free path, we obtain a condition on the diffusion coef- 
ficient: 

HV,_L656k^r 1 

A„„_, ~ z — ^ (li) 



2 yllncuf-MxMm) ^^^^ 
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For the big particles, the Stokes regime will become domi- 
nant thus we neglect the second term and the condition needed 
is: 
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As the relevant range of particle size span several order of 
magnitudes, a good approximation of K^^.^ can be obtained by 
taking the sum of these two coefficient: 
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Again we note that these limits for the diffusion coefficient 
are independent of the planet considered. Given a global solar 
abundance for TiO in the planet, Spiegel et al. (2009) stated 
that an abundance of half the solar composition at Imbar would 
be necessary to produce an observable stratosphere. Applying 
formula 26 with xum = 0.5 and Pum = 1 mbar and assuming 
a well mixed layer below 1 bar (i.e. Pi = I bar) we obtain 
Kzz,,^ = 1.4xl0'*m2s"' fora = 0.1 ^m, = 1.4xl05m2s"' 
for a = 1 |im and K^zj^^ = 1.6 X lO^m^s"' for a = 10 \iia. 
These results are of the same order of magnitude as the ones 
found by Spiegel et al. (2009). This is expected, since we com- 
pare two similar mechanisms: settling and diffusion of particles. 
However, this similarity of the results shows that the day-night 



cold trap (not considered by Spiegel et al.) is at least as impor- 
tant as the vertical one in hot- Jupiter atmospheres. Moreover, 
the condition on K^z^i^ we derived are independent of the planet 
studied and holds for very hot Jupiters such as WASP- 12b or 
WASP-33b where the vertical cold trap would be inefficient or 
nonexistent. 



5. Effective vertical diffusion coefficient 

As described previously, ID models have been extensively used 
to investigate chemistry and vertical structure of hot Jupiter at- 
mospheres, with the 3D dynamics parameterized as vertical eddy 
diffusion with a specified diffusivity (e.g Spiegel et al. 2009; 
Zahnle et al. 2009b,a; Youdin & Mitchell 2010; Line et al. 2010, 
201 1; Madhusudhan et al. 201 1; Moses et al. 201 1). 

There is no theoretical reason for the mixing by the global 
flow patterns in hot Jupiter atmospheres to behave like a one- 
dimensional diffusion process. However, deriving an a posteriori 
effective diffusion coefficient that describes as closely as possi- 
ble the averaged vertical mixing within the atmosphere can be 
a useful way to roughly characterize the strength of the verti- 
cal fluxes of material and guide ID modelers in their choice of 
vertical mixing parameters. 

A first way to define a vertical mixing coefficient from our 
simulation is to choose the that best reproduces the planet 
averaged tracer profiles. The ID model developed in Sect. 4 de- 
scribes the equilibrium between vertical diffusion of tracers with 
a specified height-dependent K^^ and their settling on the night- 
side. We tune the diffusivities to obtain a good match between 
the solutions of our ID model and the horizontal-mean tracer 
abundance versus pressure from the 3D models. To use equation 
17 we must specify the temperature of the atmosphere, constant 
in the analytical model. The temperature appears in the expres- 
sion for the Knudsen number and in the expression for the vis- 
cosity of hydrogen. Both quantities are related to the settling of 
the particles. Thus the temperature to consider is the nightside 
temperature. In our GCM nightside temepratures range from 600 
to s; 1500 K and we decided to use a mean temperature of 1000 
K. However, we note that the derived value of Kzz does not de- 
pend strongly on this choice. Then using a Kzz value proportional 
to the inverse square root of the pressure (a = 0.5 in equation 
17), we obtain a remarkably good agreement between our ID 
model (dotted lines in fig. 4) and the horizontal average of the 
3D model (solid lines of fig. 4). The resulting value for the verti- 
cal mixing coefficient that best fits the different tracer field used 
in the simulation is : 



^ 5 X 10^ 2 -1 
^^zz = , m^s , 



(27) 



which is valid over a pressure range from ~1 bar to a few //bar. 
This profile is plotted in fig. 13 (red Une). Interestingly, the same 
profile of fits the 3D simulations reasonably well over a wide 
range of particle size (see fig. 4), suggesting that does not 
depend strongly on particle size for the models considered here. 
The 3D models adopt the stellar insolation and other properties 
for HD 209458b, so the results are most germane to that planet; 
the mixing rates are likely to be larger for hotter planets and 
smaUer for cooler planets than implied by Equation 27. 

Another way to define a one dimensional vertical mixing co- 
efficient from the three dimensional simulation is to find the K^^ 
that leads to an upward diffusive flux of material that matches 
the averaged vertical flux produced by the dynamics. This can 



14 



Parmentier et al.: 3D mixing in hot Jupiter atmospheres I 



be written (e.g. Chamberlain & Hunten 1987, p. 90): 



where the brackets represent the horizontal average along iso- 
bars over the whole planet. This expression does not necessitate 
any assumption on the functional form of K^z nor on the night- 
side mean temperature. It also has the advantage that no com- 
parisons or fits to a ID diffusion model are necessary; the effec- 
tive values of Kzz can be derived directly from the 3D GCM 
data via Equation 28. As a trade-off, this expression depends 
on the vertical tracer gradients that are affected by long term 
temporal variability (see Sect. 3.4), which are not smoothed out 
completely given the limited integration time of the simulation. 
This leads to some strong vertical variation of K^z- The profile of 
calculated from Equation 28 are shown in the black curves 
of fig. ??. Despite the vertical fluctuations, the overall shape of 
Kzz obtained with this method is close to the estimate using the 
planet averaged tracer profile (see the red curve in fig. 13). As 
before, we note that the derived value does not depend strongly 
on the particle size, consistent with the fact that Equation 27 al- 
lows a good fit between ID models and our 3D model for a wide 
range of particle sizes. 




Kzz (m^/s) 



Fig. 13: Vertical diffusion coefficient from the 3D model (black 
lines), from the ID fit to the 3D model (red) and the root mean 
square of the vertical velocity times the vertical scale height, a 
common estimate of Kzz in the literature (blue). 

Several previous studies have attempted to estimate the ver- 
tical diffusion coefficient in hot Jupiter atmospheres (Cooper & 
Showman 2006; Moses et al. 2011; Heng et al. 2011b; Lewis 
et al. 2010). Cooper & Showman (2006) adopted an estimate 
for Kzz based on the product of a root-mean-square vertical ve- 
locity from their 3D GCMs and an appropriately chosen ver- 
tical length scale following the formulation of Smith (1998). 
Moses et al. (2011) and Lewis et al. (2010) followed a similar 
procedure but adopted an atmospheric scale height for the verti- 
cal length scale. These estimates are crude, although Cooper & 



Showman (2006) showed that this formulation for Kzz allows ID 
models to match the full tracer profiles from 3D GCMs reason- 
ably well. More recently, Heng et al. (201 lb) used the magnitude 
of the Eulerian mean streamfunction as a proxy for the strength 
of the vertical motions and derived a vertical mixing coefficient 
of Kzz ~ lO^m^s"'. Note, however, that the Eulerian-mean ve- 
locities are known to be a poor descriptor of tracer advection 
rates in planetary atmospheres, since eddy mixing often domi- 
nates over transport due to the Eulerian-mean circulation (see, 
e.g. Andrews et al. 1987, Chapter 9). 

Although we confirm that mixing in hot- Jupiter atmospheres 
is strong, we find a value that is significantly smaller than the 
previous ones. In particular, our value is two order of magni- 
tude smaller than what is obtained when multiplying the vertical 
scale height by the root mean square of the vertical velocity (blue 
curve in fig. 13), a common estimate for Kzz in the literature. 

As stated in Sect. 2.1, the model does not include any sub- 
grid vertical diffusion coefficient. Yet, given the huge values for 
Kzz that we derive from the resolved flow, it seems unlikely that 
sub-grid turbulent mixing would contribute significantly to the 
global mixing. However, the interaction between small-scale tur- 
bulence and the global flow might not be trivial and a more de- 
tailed study would be required to draw a firm conclusion. 

We emphasize that mixing by the global circulation appears 
to be planet- wide and differs from region to region. Although the 
globally averaged dynamics seem to be reasonably described by 
a vertical mixing coefficient, that is not the case for the local 
flow in the simulation. It is therefore difficult to define mixing 
coefficient values for particular locations in the planet. 

Along the equator, where the strong flow efficiently mixes 
the tracers longitudinally, we expect a good agreement of our 
ID model to the 3D flow. Indeed, using the value of Kzz de- 
rived in Sect. 5 we realize that the spread of the tracer profiles 
along the equator (fig. 14) is of the same order of magnitude 
as the spread predicted by the ID model (middle panel of fig. 
12). However, in the ID model, profiles from equally sampled 
longitude are equally spaced in abundances whereas the profiles 
obtained from the 3D model are sometimes packed together and 
sometimes widely spread, denoting an unequal strength of the 
vertical mixing longitudinally. 

6. Applications 

6. 1 . Presence of a stratosphere on hot Jupiters 

The presence of TiO in hot- Jupiter atmospheres is a leading hy- 
pothesis for the absorber needed to create temperature inver- 
sions (Hubeny et al. 2003; Fortney et al. 2008). However, TiO 
condenses at temperatures lower than » 2000 K. On most hot 
Jupiters, while gaseous TiO can be stable on the dayside, it 
should condense on the night side. Thus, our results allow us 
to address the question of whether TiO can remain suspended in 
the atmospheres of hot Jupiters, and hence whether TiO-induced 
stratospheres are indeed viable. Spiegel et al. (2009) predicted 
that for a solar abundance of TiO in the planet, an abundance of 
0.5 times the deep abundance at 1 mbar is necessary to main- 
tain a temperature inversion on the day-side. Our simulations 
suggest that if TiO condenses into particles bigger than several 
micrometers, the day-night cold trap will be sufficiently efficient 
to deplete it from the dayside. If, on the contrary, TiO cannot 
condense in particles bigger than several micrometers, it should 
remain present on the dayside and produce a stratosphere. 

Spiegel et al. (2009) emphasized the deep cold trap in the 
deep stratified layer at pressures exceeding tens of bars where 
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Tracer abundance 

Fig. 14: Time averaged tracer abundance along the equator. The 
profiles are equally spaced in longitude. The dayside profiles are 
in red whereas the nightside profiles are in blue. 



dynamical mixing rates are probably low. In planets that exhibit 
such a cold trap, mixing TiO upward to altitudes where it could 
be affected by the day-night cold trap may be difficult. On the 
other hand, the opacities and therefore temperature structure in 
these deep regions are rather uncertain; moreover, hot Jupiters 
that are particularly highly irradiated would exhibit warmer tem- 
peratures and would therefore be less likely to exhibit such a ver- 
tical cold trap. Some very hot Jupiters may therefore lack such a 
deep cold trap entirely. On these planets, the day-night cold trap 
at low pressure (P < 1 bar) would then dominate. 

The spatial variations in the gaseous tracer abundance in our 
3D models (fig. 5) suggest that the TiO abundance on the day- 
side, and hence the stratosphere itself, could be patchy, with 
some regions of the dayside exhibiting a stronger temperature 
inversion than others. This would have interesting consequences 
for the interpretation of dayside infrared spectra. 

The time variability described in Sect. 3.4 could affect the 
presence of the stratosphere, leading to strong temporal vari- 
ability in the upper atmospheric temperatures. In particular, the 
tracer abundance averaged over much of the dayside exhibits 
large-amplitude fluctuations (fig. 8), particularly for larger par- 
ticle sizes. This suggests that, at least for some planets, the TiO 
abundance could fluctuate between values large enough to gen- 
erate a stratosphere and values too small for a stratosphere to 
form. The stratosphere itself might then fluctuate episodically 
in time. Although not included in our current models, there is 
the possibility of feedbacks with the flow itself, since the pres- 
ence (or absence) of a stratosphere exerts a significant impact 
on the flow structure and vertical mixing rates. If the feedback 
is positive, i.e. if the presence of high temperatures in the upper 
atmosphere enhances the mixing, then hot Jupiters could oscil- 
late between a state with strong vertical mixing and stratospheric 
heating by TiO and a state with no stratospheric heating and less 
vertical mixing. This is a two-state atmosphere analogous to that 
described by Hubeny et al. (2003). However, this possibility re- 



mains speculative and further models that include the feedback 
of the tracer field on the flow are necessary to draw a firm con- 
clusion. 

TiO is thought to be the major Ti-bearing gas in hot- Jupiter 
atmospheres (see Lodders 2002). However, other Ti-bearing 
gases, Ti02 being the most abundant, are believed to be the 
condensable Ti-bearing species. Thus, a parcel of gas experi- 
encing a sudden drop in temperature due to its advection to the 
night side might not see all its titanium incorporated into con- 
densate, but rather only the titanium atoms that already reside 
in other Ti-bearing molecules such as Ti02(g). The reaction be- 
tween gaseous species TiO -i- H2O ^ TiOa + H2 is fast under the 
conditions relevant to the dayside of hot Jupiters (see Fortney 
et al. 2009). However, the reaction might be kinetically inhib- 
ited on the night side of the planet where the temperature drops 
significantly, leading to a smaller depletion of Ti than in the hy- 
pothetical case where TiO(g) could condense by itself. 

The dayside abundance of tracers and thus any conclusion on 
the possibility to maintain a stratosphere in hot Jupiter depends 
on the size of the condensate, a free parameter in our study. This 
size results from complex microphysical processes. Once on the 
night side, TiO is over-saturated, and thus Ti-bearing conden- 
sates are expected to appear We used the formalism of Woitke & 
Helling (2003) to calculate the characteristic growth time scale 
of TiOa particles assuming that all the titanium is contained in 
Ti02(g) and that this last fully saturates the atmosphere. We ex- 
pect that condensate will hardly form if the condensate growth 
time scale exceeds the advective time scale, i.e., the time for the 
jet to travel across one hemisphere. Thus all the particles above 
the black line in fig. 15 are unlikely to form. The low elemental 
abundance of titanium — solar abundance is 10"^ compared to H 
(Lodders 2002) — ^kinetically inhibits the formation of microme- 
ter size particles above lOmbar. 

However, titanium is not the only element that can form con- 
densates on the night side of hot Jupiters. Silicates are believed to 
condense and could incorporate titanium atoms into their grains. 
Sub-micron size Ti02 particles could even be used as seeds for 
the formation of silicates grains. In that case the relevant time 
scale is not tied to the growth of Ti02 grains but rather to the 
growth of Si02 based grains (MgSi02 for example). We calcu- 
lated the growth time scale of Si02-based grains using the same 
formalism as for Ti02 grains, assuming that Si02(g) is fully sat- 
urated and that all the silica atoms are in Si02(g) molecules. As 
can be seen in fig. 16, particles as big as 10 \xm can form at pres- 
sures as low as lOmbar in that case. 

The estimates in fig. 16 suggest that the growth time of 
siUcate particles is comparable to advection times across a 
hemisphere in the pressure range where the stratosphere forms 
(~0.01-1 mbar). This comparison suggests that more detailed 
calculations, coupling the 3D dynamics to microphysics that al- 
low self-consistent prediction of particle growth, may be neces- 
sary to obtain a firm conclusion about whether particle growth 
timescales are sufficiently long to inhibit loss of TiO from the 
atmosphere. 

6.2. Clouds in hot Jupiter atmosplieres 

The huge day-night temperature difference and cold nightside 
temperatures predicted on many hot Jupiters at low pressure 
(e.g. Showman et al. 2009) suggest that, in addition to TiO, a 
wide range of other chemical species, including silicates and 
iron, will condense on the nightside. Some of them could also 
stay in a condensed state in part or all of the dayside hemi- 
sphere. The Rayleigh scattering slope in the transmission spec- 
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Fig. 15: Approximate growth time scale of TiOi grains in HD 
209458b. The bottom black line shows the advective timescale 
and the top one is ten times the advective timescale. 
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Fig. 16: Approximate growth time scale of SiOa grains in HD 
209458b. The bottom black line shows the advective timescale 
and the top one is ten times the advective timescale. 



trum of HD 189733b, first observed by Lecavelier Des Etangs 
et al. (2008) and later confirmed by numerous observations (see 
Pont et al. (2012) for a review of the different observations of this 
planet) is best fitted by models including sub-micron sized par- 
ticles. Then the strong spatial and temporal variations observed 
in our model can also be interpreted as spatial variation of the 
cloud coverage in the atmosphere of hot- Jupiter This could lead 
to albedo variations along the dayside of the planet. Due to the 
eastward displacement of the hottest point of the planet, clouds 
might be more able to form west of the substellar point than 
east of it, leading to a longitudinal variation of the albedo that 
contributes to the spatial variation described before. The large 
amount of data in the visible from the Kepler space telescope 
is ideal to search for such a spatial and temporal variability in 
albedo pattern of tidally locked planets. 



6.3. Parameter retrieval 

Atmospheric characteristics of hot Jupiters are usually de- 
rived from disc-integrated fluxes (for secondary eclipses mea- 
surements) or limb-integrated transmission (for transits spec- 
troscopy). Therefore, interpretation of the data is usually done 
using one dimensional atmospheric models, assuming an ho- 
mogenous atmosphere, both in term of temperature structure and 
composition. However, given the strong spatial variability in the 
tracer distribution both on the dayside hemisphere and in the 
limb profiles observed, some exoplanet spectra might be better 
understood by considering spatial variation in the atmospheric 
profiles and chemical composition along the planet. For exam- 
ple, an inhomogeneous distribution of TiO would lead to strong 
brightness differences in the emitted flux from the different lo- 
cations of the planet, which could affect features such as the ap- 
parent eastward offset of the brightest spot of the planet. Hazes 
could also show a similar behavior, leading to planetary spec- 
tra that might be better explained with the combination of two 
different (cloudy and cloudless) ID models — as is the case in 
recent models of brown dwarfs (e.g. Marley et al. 2010). The 
rising technique of phase mapping (see Majeau et al. 2012; De 
Wit et al. 2012) might soon allow us to constrain better the spa- 
tial inhomogeneities in the disk of the planet. 

Conclusion 

We presented global, three-dimensional numerical simulations 
of the atmospheric circulation of HD 209458b including a pas- 
sive tracer that can be either gaseous or solid. This is the 
first circulation model of a hot Jupiter to include a representa- 
tion of clouds. We applied our model to chemical species that 
are gaseous on the dayside but condense on the nightside of 
the planet. Given the strong day/night contrast present in hot 
Jupiters, our model applies to a wealth of different chemical 
species such as titanium and vanadium oxides and silicate ox- 
ides among others. 

We demonstrate that, although hot Jupiter atmospheres are 
believed to be stably stratified (i.e. locally non-convective), they 
are strongly mixed. In the presence of a background gradient 
of chemical species, large-scale circulation patterns naturally 
create upward mixing. This mixing, resolved by the GCM, is 
strong and likely dominates over molecular, convective or tur- 
bulent mixing. In HD 209458b, the mixing is strong enough to 
keep a condensable species lofted if it condenses into particles 
smaller than a few microns on the night side of the planet. 

The coupling between 3D flow and particle settling leads 
to strong spatial and temporal variations in the abundance of a 
given chemical species. Around 0.1 mbar, the tracer abundance 
is homogeneous in longitude but exhibits a large latitudinal vari- 
ation, the equator being more depleted than the poles. Around 
1 mbar, at high latitudes, the day/night contrast becomes im- 
portant. According to our models, variability of up to ~50% in 
the dayside tracer abundance, and of up to ~75% in the tracer 
abundance along the limb, can occur for sufficiently large parti- 
cle sizes (~5 /im). This variability characteristic periods ranging 
from days to ~50-100 days. 

These results can be applied to a wide range of molecules in 
hot Jupiters atmospheres. Titanium oxide, the best candidate for 
creating a temperature inversion in the dayside of hot Jupiters, 
should condense on the night side of most planets. Our results 
imply that the day/night cold trap could impede the formation 
of a stratosphere in the dayside if TiO condenses into particles 
bigger than a few microns on the night side. Growing particles 
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to such a size seems difficult when TiO alone is considered due 
to its small abundance. However, TiO can be incorporated into 
condensates from more abundant gases such as sihcate oxides. 
In that case the day/night cold trap could be strong enough to im- 
pede the formation of a hot stratosphere on the day side. Spatial 
variabiUty of TiO could significantly affect the dayside temper- 
ature structure and exert interesting efi'ects on infrared spectra 
and lightcurves. For example, mid-to-high latitudes might keep 
enough TiO to create an inversion whereas the equator, more de- 
pleted, might not be able to sustain the inversion. Such a latitu- 
dinal contrast could be observed using the phase mapping tech- 
nique (e.g. Majeau et al. 2012; De Wit et al. 2012). The temporal 
variability observed in the model could lead to the appearance 
and disappearance of the stratosphere on timescales of ~10- 
100 earth days. Highly irradiated planets have significant ther- 
mal emission in the Kepler bandpass (e.g. Spiegel & Burrows 
2010, for Hat-P-7b). Using the long photometric series from the 
Kepler Space Telescope, such a variability might be observable. 

Our results also apply to silicate hazes. Temporal and spatial 
variabihty in the cloud coverage could strongly aff'ect the albedo 
and the thermal emission of the planet. For moderately irradiated 
planets, the Kepler spacecraft observes the reflected light of the 
star by the planet (e.g. Demory et al. 201 1, for Kepler-7b). Thus, 
the time series from Kepler could be used to build albedo, and 
therefore cloud maps of the planet. 

Although there is no theoretical reason for the upward mix- 
ing driven by the global circulation to be diffusive, it is inter- 
esting to quantify the averaged vertical mixing with a diffusive 



model. The value K„ value of Kr, 
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m^s , between 



~1 bar and ~1 mbar, can be used in ID models of HD 209458b. 

This study, the first one including a representation of clouds 
in a GCM of hot- Jupiter, shows that hot- Jupiter atmospheres are 
strongly mixed and that large scale spatial and temporal vari- 
ability are expected in any condensable chemical constituents. 
Today, observers can already detect longitudinal variations in the 
emitted thermal flux of the planet. In the next decade, both longi- 
tudinal and latitudinal variations in thermal emission and albedo 
of the hot Jupiters will be observable, expanding the study of 
weather to extra-solar planets. 
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Appendix A: Departure from the Cunningham 
velocity 

The Stokes-Curmingham velocity defined in equation 3 is de- 
rived under the assumption of low Reynolds number. Therefore 
it is not valid for turbulent flow and other expressions may be 
used when the Reynolds number increases. Here we derive bet- 
ter laws for intermediate and large Reynolds number. 

A.1. Low Knudsen number 

For small Reynolds number and small Knudsen number, the drag 
force exerted by a fluid on a sphere at rest is considered propor- 
tional to the kinetic energy of the fluid and the projected area of 
the sphere. The coefficient of proportionality, or drag coefficient, 
Cd is given by : 

Co = (A.1) 

Then, equating gravity and drag forces leads to the settling ve- 
locity of a particle in an atmosphere : 



(A.2) 



Where pp is the density of the particle. For smaU Reynolds 
numbers and high Knudsen number, Co = 24 is constant and 
the settling velocity is the Stokes velocity. When increasing the 
Reynolds number, the non linear terms of the Navier-Stokes 
equation become important and Co is no longer constant. We 
used tabulated values of the drag coefficient as a function of the 
Reynolds number given by Pruppacher & Klett (1978). We as- 
sume that Cd = 24 when the Reynolds number reaches 1 to stay 
consistent with Stokes ffow and that Cd reaches its asymptotic 
value, Cd = 0.45, when Nr^ = IR^ = 1000 and fft the relation- 
ship : 

logioCA^R.) =-1.215047 + 0.923242 logioCCoAf^) 

-0.031293 logio(Cz,<^)2 (A.3) 

Then we foUow the same method as Ackerman & Marley 
(2001). Noting that : 



2 ^2pga^(pp-p) 
CdN^ = J^f 



(A.4) 



is independent of the velocity, we use the fit of equation A.3 and 
extract the velocity : 



^ JLlO-l-^'f^^^rf^a^V"'"'' (A.5) 



2pa \ 'irf- j \ 2>Tf 



A.2. High Knudsen number 

In the free-molecular regime, calculations have been made by 
Probstein (1968) leading to an expression for the drag coeffi- 
cient: 



2 _ 2s^ + 1 , , ' 
Cd = ^ + exp(-sl) + ■ 



24 



erfis,) (A.6) 



where Sa is the ratio of the object velocity over the thermal speed 

of the gas (Vt = yj^^ - "T^ where c is the sound speed) and 
erf is the error function. 

For velocities much smaller than the sound speed, Sa — > 
and we can use an equivalent of the error function in 0: 



2 2 2 
erfis^) = ——e-Hsa + r^a) + ^(4) 



(A.7) 



Using this equation inside Equation A.6 and taking the hmit 
Sa 0, the term e"^" goes to 1 and the terms proportional to ^ 
cancels out leading to: 



16 \ Vt 
3^/ V 



3 3 V?r/ 

For velocities much greater than the sound speed, the limit 
of equation A.6 when — > 0° is: 

Cz, ~ 2 (A.9) 

In order to simplify equation A.6 we use the following ex- 
pression for the drag coefficient at high Knudsen number : 



3 -svi/v- <^-"'> 

Our approximation ffts correctly the exact expression in the 
limit of low and high velocities. In between the difference to the 
exact expression is at most 30%. Replacing Cd by its value in 
equation A.2 we obtain a second order equation for the velocity 

^ / y \^ /2V7r 16 \ y 8 Ap „ 



which leads to : 



/, 96 Vstokes „ 
1 1 H ; — . —TZ r^N ■ 



with A = + 3^)' ^^^^ th^ speed becomes smaU com- 
pared to the sound speed (ygtokes *^ ^t) we obtain : 
96 



Vf 



8 V^A 



^N^^stokes ~ l-61^rNystokes 



(A.13) 



which is in good agreement with equation 3, derived for high 
Knudsen number and small Reynolds numbers. 

A.3. Comparison with the Cunningham velocity 

Figure A.l shows the ratio of the Cunningham velocity (see 
equation 3) to the ones we just derived. The difference is no- 
ticeable only for particles of the order of 100 |am at pressures 
less than the 10"'' bar level and exceeding the 10 bar level. This 
difference is always less than one order of magnitude and con- 
cern only a tiny portion of the parameter space which has little 
relevance to our study (the largest particle sizes considered in 
our 3D models is lOyum). Thus we decided to neglect this dis- 
crepancy in the main study. 
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Fig. A.l: Ratio of the Cunningham velocity to the more sophis- 
ticated model of the particle settling velocity considered in the 
Appendix. 
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